clear

%set a name for the estimation output
estimation_name='abc_rn';

%set risk neutral or not
risk_neutral=1;         %if set to zero, takes the individual lottery choices for crra utility
                        %if set to one, assumes risk neutrality
                        %if set to two, estimates r as a free parameter
                        %if set to three, assumes log-utility (r=1)
                        
%set constraints
alpha_zero=0;           %sets alpha to zero
beta_zero=0;            %sets beta to zero
kappa_zero=0;           %sets kappa to zero
delta_zero=1;           %sets delta to zero
gamma_zero=1;           %sets gamma to zero

alpha_negative=0;       %choose -1 to set alpha positive, 1 to set alpha negative, 0 for no constraint on the sign
beta_negative=0;        %choose -1 to set beta positive, 1 to set beta negative, 0 for no constraint on the sign
kappa_negative=0;       %choose -1 to set kappa positive, 1 to set kappa negative, 0 for no constraint on the sign
delta_negative=0;       %choose -1 to set delta positive, 1 to set delta negative, 0 for no constraint on the sign
gamma_negative=0;       %choose -1 to set gamma positive, 1 to set gamma negative, 0 for no constraint on the sign
lambda_negative=-1;     %choose -1 to set lambda positive, 1 to set lambda negative, 0 for no constraint on the sign


%import choices from experiment
filename='mnl_data.csv';
dt=csvread(filename);
file2='risk.csv';       % contains the Eckel-Grossman lottery choices
risk=csvread(file2);

no_trials=18;
n = length(dt)/no_trials; % n is the number of subjects in the data

options = optimoptions(@fmincon,'MaxFunEvals',6000);
n_steps = 6;           % number of starting points used per parameter

Results = [];
AllResults = [];

if risk_neutral==2
    n_params = 7;
    r_steps = n_steps;
else
    n_params = 6;
    r_steps = 1;
end

%sets with initial values
if alpha_negative == -1
    a_min = 0;
    a_max = 1;
elseif alpha_negative == 1
    a_min = -1;
    a_max = 0;
else 
    a_min = -1;
    a_max = 1;
end

if beta_negative == -1
    b_min = 0;
    b_max = 1;
elseif beta_negative == 1
    b_min = -1;
    b_max = 0;
else 
    b_min = -1;
    b_max = 1;
end

if kappa_negative == -1
    c_min = 0;
    c_max = 1;
elseif kappa_negative == 1
    c_min = -1;
    c_max = 0;
else 
    c_min = -1;
    c_max = 1;
end

if delta_negative == -1
    d_min = 0;
    d_max = 1;
elseif delta_negative == 1
    d_min = -1;
    d_max = 0;
else 
    d_min = -1;
    d_max = 1;
end

if gamma_negative == -1
    e_min = 0;
    e_max = 1;
elseif gamma_negative == 1
    e_min = -1;
    e_max = 0;
else 
    e_min = -1;
    e_max = 1;
end

if lambda_negative == -1
    l_min = 0.01;
    l_max = 14.01;
elseif lambda_negative == 1
    l_min = -14.01;
    l_max = -0.01;
else 
    l_min = 0.01;
    l_max = 14.01;
end

r_min = -0.1; r_max = 1.6;

A_unc = zeros(n_params,n_params);
b_unc = zeros(1,n_params);
Aeq_unc = zeros(n_params,n_params);
beq_unc = zeros(1,n_params);

if alpha_zero == 1
    Aeq_unc(1,1) = 1;
    a_steps = 1;
else
    a_steps = n_steps;
end
    
if beta_zero == 1
    Aeq_unc(2,2) = 1;
    b_steps = 1;
else
    b_steps = n_steps;
end

if kappa_zero == 1
    Aeq_unc(3,3) = 1;
    c_steps = 1;
else
    c_steps = n_steps;
end

if delta_zero == 1
    Aeq_unc(4,4) = 1;
    d_steps = 1;
else
    d_steps = n_steps;
end

if gamma_zero == 1
    Aeq_unc(5,5) = 1;
    e_steps = 1;
else
    e_steps = n_steps;
end

A_unc(1,1) = alpha_negative;
A_unc(2,2) = beta_negative;
A_unc(3,3) = kappa_negative;
A_unc(4,4) = delta_negative;
A_unc(5,5) = gamma_negative;
A_unc(6,6) = lambda_negative;

risky2=[4.71 2.95 1.19 0.77 0.32 -0.13];    % r parameters from lotteries in the first session
risky=[1.61 1 0.39 0.25 0.08 -0.09];        % r parameters from lotteries in all other sessions

ind_data = zeros(no_trials,11,n);

parfor i = 1:n

    m1 = no_trials*(i-1)+1; m2 = no_trials*i; id = dt(m1,1);
      
    subdt = [dt(m1:m2,3:6),dt(m1:m2,7:9),dt(m1:m2,10),dt(m1:m2,11),risk(i)*ones(no_trials,1),id*ones(no_trials,1)];
    ind_data(:,:,i)=subdt;
    ids(i)=id;
end

parfor i = 1:n
     
    f_unc=@(pm2)lolik_indiv(pm2,ind_data(:,:,i),risk_neutral);
    
    ind_Results = [];
    for a=1:a_steps
        for b=1:b_steps
            for c=1:c_steps
                for d=1:d_steps
                    for e=1:e_steps
                        for l=1:n_steps
                            for r=1:r_steps
                        
                            try
                             
    a_start=a_min+a*((a_max-a_min)/(a_steps+1));
    b_start=b_min+b*((b_max-b_min)/(b_steps+1));
    c_start=c_min+c*((c_max-c_min)/(c_steps+1));
    d_start=d_min+d*((d_max-d_min)/(d_steps+1));
    e_start=e_min+e*((e_max-e_min)/(e_steps+1));
    l_start=l_min+l*((l_max-l_min)/(n_steps+1));
    
    if risk_neutral == 2
        r_start=r_min+r*((r_max-r_min)/(r_steps+1));
        pm0 = [a_start b_start c_start d_start e_start l_start r_start];
    else
        pm0 = [a_start b_start c_start d_start e_start l_start];
    end
    
    [pm2,fval] = fmincon(f_unc,pm0,A_unc,b_unc,Aeq_unc,beq_unc,[],[],[],options);

    alpha_est = pm2(1); beta_est = pm2(2); kappa_est = pm2(3); delta_est=pm2(4); gamma_est=pm2(5); lambda_est=pm2(6);

    if risk_neutral == 2
        rho_est = pm2(7);
        rho0 = pm0(7);
    elseif risk_neutral == 0
        if ids(i)<200
            rho_est = risky2(risk(i));
            rho0 = risky2(risk(i));
        end
        if ids(i)>=200
            rho_est = risky(risk(i));
            rho0 = risky(risk(i));
        end
    elseif risk_neutral==3
        rho_est = 1;
        rho0 = 1;
    else
        rho_est = 0;
        rho0 = 0;
    end
    
    ind_Results = [ind_Results;ids(i), alpha_est, beta_est, kappa_est, delta_est, gamma_est, rho_est, lambda_est, -fval, pm0(1), pm0(2), pm0(3), pm0(4), pm0(5), rho0, pm0(6)];

           catch ME
    
    ind_Results = [ind_Results;ids(i), 9999, 9999, 9999, 9999, 9999, 9999, 9999, -9999999999, 9999, 9999, 9999, 9999, 9999, 9999, 9999];
    
            end

                            end
                        end
                    end
                end
            end
        end
    end
    
    logl=ind_Results(:,9);
    [M,I]=max(logl);
    
    Results = [Results;ind_Results(I,:)];

end    

%save workspace
save(fullfile('output','indiv','mat_files',estimation_name));

% export results
id_est = Results(:,1);
a_est = Results(:,2);
b_est = Results(:,3);
c_est = Results(:,4);
d_est = Results(:,5);
e_est = Results(:,6);
r_est = Results(:,7);
l_est = Results(:,8);
ll_est = Results(:,9);   
a0_est = Results(:,10);
b0_est = Results(:,11);
c0_est = Results(:,12);
d0_est = Results(:,13);
e0_est = Results(:,14);
r0_est = Results(:,15);
l0_est = Results(:,16);
T = table(id_est,a_est,b_est,c_est,d_est,e_est,r_est,l_est,ll_est,a0_est,b0_est,c0_est,d0_est,e0_est,r0_est,l0_est);
writetable(T,fullfile('output','indiv','estimates',estimation_name));

